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The  Equation  of  State  for  Ammonia 

by 

Lester  Haar  and  John  Gallagher 

1.  Introduction. 

In  this  report,  we  present  an  outline  of  the  basic  results  of  an 
extensive  correlation  of  the  thermodynamic  properties  of  ammonia^.  The 
purpose  of  this  report  is  to  present  an  early  working  version  of  that  correlation 
for  the  convenience  of  those  workers  with  an  immediate  need 
for  such  results  and     to  which  they  can  refer  in  an  unambiguous  manner. 
This  report  has  been  prepared  in  answer  to  several  urgent  requests  for 
such  a  version. 

We  present  the  basic  equations  used  for  the  description  of  the 
properties  of  ammonia  and  listings  of  computer  programs  from  which  all  the 
thermodynamic  properties  for  ammonia  can  be  calculated  for  the  temperature  range 
including  the  triple  point  temperature  to  about  5/3  the  critical  temperature  and  for 
the  pressure  range  including  the  dilute  gas  to  about  8000  bars.  The 
reference  state  for  all  properties  is  the  ideal  gas  at  zero  kelvin.  The 

physical  constants  used  are  consistent  with  those  recommended  by  Cohen  and 
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Taylor  .     The  mass  of  a  mole  of  ammonia  was  taken  to  be  17.03026  grams  . 

2.  The  derivation  of  the  thermodynamic  surface. 

In  this  section,  we  describe  the  results  of  a  fit  of  selected  P,p,T 
data  schematically  represented  in  figure  1  to  an  analytic  equation  of  state  in 
the  least  squares  sense.     A  detailed  discussion  of  the  process  of  data  selection 
and  of  the  selection  of  parts  of  the  analytic  equation  is  contained  in 
reference  1.     The  equation  so  obtained  contains  the  pressure  as  a  44  term  double 
power  series  function  of  temperature  and  density.     This  equation  can  then  be  used  to 


reproduce  all  the  available  P,p,T  experimental  data  as  well  as  to  produce 
limited  extrapolations  (based  on  thermodynamic  arguments)  of  the  surface  into 
important  regions  where  data  are  sparse.     The  range  of  the  equation  is  bounded 
at  low  temperatures  by  the  triple  point  temperature  (195. 48K)  and  the  melting 
curve  for  the  liquid,   and  at  high  temperatures  by  the  isotherm  at  750K  (which 
is  approximately  5/3  the  critical  temperature).     The  pressure  range  extends 
to     8,000  bar.     No  mathematical  constraints  were  imposed  on  the  equation ,  and 
only  P , p , T  data  were  used  in  the  least  squares  f it . 

Following  Keen  an  et  al.     the  Ilelmholtz  free  energy  function  was 
represented  in  reference  1  as  the  sum  of  two  terms:   the  first  is  the 
contribution  from  the  equation  of  state;   the  second  is  a  function  of 
temperature  only  and  refers  to  the  properties  of  the  ideal  gas.     Thus  the 
Helmholtz  free  energy  was  expressed, 

A(p,T)  «  A(p,T)  +  A°(T),  (D 

where  A°(T)   is  the  contribution  of  the  ideal  gas.     T  is  the  absolute 
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temperature  in  Kelvin  and  p  is  the  density  in  grams  per  cm  .       A  quantity 
0(p,T)  was  defined  by 

A(p,T)  =  RT  [In  p  +  pQ  (p,T)]  (2) 

■  ■ 

Since  P  -  p2  3A/dp, 

Eqs.  (l)and(2)  yield, 


P  -  pRT  [1  +  pQ  +  p2  3Q/9p]  . 


(3) 


We  note  that 


Q  (p-0)  -  B2,  (3a) 


where        is  the  second  virial  coefficient.     The  form  chosen  for  Q  was 
9  6 

Q  "     7         7  P1"1  (t-t J5'1  ,  (4) 


i-1  j-1 


where  x  =         ,        =  1.2333498  and  R  is  the  gas  constant.     Equation  3  for  the 
pressure  was  fitted  in  the  least  squares  sense  to  the  experimental  P,p5T  data, 

The  results  of  this  fit  are  values  for  the  constants  a.,   listed  in  table  1. 

<  t     <•  ij 

(a..  =  0  for  all   ..  not  listed  in  table  1).     By  differentiation  of  eq  (1) 
we  obtain  for  the  entropy, 


S(p,T)  -  R[Jtn  p  +  pQ  -  pt  9Q/9t]  +  S°(T),  (5) 

the  internal  energy, 

E(p,T)  -=  RpTt3Q/8t  +  E°(T),  (6) 


the  constant  volume  heat  capacity, 


C  (p,T)  -  -  Rpr2  32Q/8T2  +  C  °, 
v  v 


(7) 
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Table  1 


i 

j 

a,. 

* 

1 

1 

1 

2 

-1j.z95oz jo/j 

1 

** 

3 

-8. 121177091!) 

1 

A 

-6.9690043553 

1 

5 

-9.73658023A9 

1 

6 

3. A8166A2617 

2 

1 

8.8100AA5762 

2 

2 

-5.07895A8707 

2 

3 

-68.261583422 

2 

A 

-74. 727156949 

2 

5 

49.751854179 

2 

6 

-14.487156374 

3 

1 

-10.467902857 

3 

2 

361.91907645 

3 

3 

1327. 8270222 

3 

A 

1484.2843304 

3 

5 

-82.229122939 

3 

6 

20.170856719 

A 

1 

75.049574001 

4 

2 

-2103. 9451938 

A 

3 

-7576.1007937 

A 

A 

-8334.8746422 

A 

5 

43.998475959 

A 

6 

-9.2773376718 

5 

1 

-409.02964153 

5 

2 

6212.2822515 

5 

3 

22341.800329 

5 

A 

23618.791735 

6 

1 

1072.479955 

6 

2 

-10816.10642 

6 

3 

-38259.344112 

6 

A 

-38233.534003 

7 

1 

-1471.4013145 

7 

2 

11195.138723 

7 

3 

38544.628190  " 

7 

A 

35887.294649 

8 

1 

1046.2341301 

8 

2 

-6365.7466698 

8 

3 

-21314.815310 

8 

A 

-18162.094974 

9 

1 

-305.80081169 

9 

2 

1532.0616045 

9 

3 

5021.6962092 

9 

A 

3812.3691534 

the  enthalpy  function, 


H(p,T)  -  RT[pQ  +  p2  3Q/3p  +  pT  3Q/3t  +  1]  +  E°  (T)  ,  (8) 
the  constant  pressure  heat  capacity, 

Cp  (p.T)  -    Cv  +  R.  |  ,  (8a) 

where, 

a  -  1  +  pQ  +  p2  3Q/3p  -    pT  3Q/3T  -  p2T  32Q/3x3p  , 

and 

3-1+2  pQ  +4  p2  3Q/3p  +  p32Q/3p2, 
the  heat  capacity  for  the  saturated  fluid  , 

T      (3P/3T)  dPe 
CS      UP      ^    (3P/3p)T  dT     *  (9) 

where  Pg  is  the  vapor  pressure  of  the  liquid  • 

S°,  E°  and  C^0  are  the  corresponding  contributions  obtained  from  A°(T) 
A°(T)  -  (G°  -  E°)  -  RT(1  +  In  A. 8180  T)  , 

*°  >  >  df  A°<«  • 


Eo.Ao(I)  .TdA°W.  (10) 
di 


C  °  -  -  T  d2  A°(T)/dT2, 
v 


where  G°  is  an  analytic  representation  of  the  results  reported  by  Haar5  for 
the  properties  of  the  ideal  gas  state, 
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G  -  E 


=  a,   in  T  a±  T1"3,  (10a) 


RT  1 

i=2 


where  E     is  the  energy  for  the  ideal  gas  at  0  K.     The  coefficients  a.  are 
o  1 

listed  in  table  2. 

Table  2 


al 

-3.872727 

a2 

.64463724 

a3 

3.2238759 

*4 

-.0021376925 

a5 

.86890833xl0-5 

a6 

an 

-.24085149xl0~7 

a7 

.36893175xl0"10 

a8 
a9 

-.35034664xl0"13 

,2056303xl0"16 

aio 

.6853420xl0_2° 

all 

.9939243x10" 24 

The  heat  capacity  values  and  the  other  thermodynamic  functions  calculated 
from  eq   (10a)   for  the  temperature  range  100  K  <  T  <  1000  K  agree  with  those 
tabulated  in  reference  5  to  within  the  accuracy  of  those  values. 


7 

Though  equations   (l-10a)  are  complete,     it  is  necessary  to  introduce  the 

Gibbs  phase  conditions  in  order  to  calculate  the  properties  for  the  coexisting 

phases.     However,  it  was  shown  by  Haar  and  Gallagher^,   that  almost  negligible 

error  results  if  an  explicit  relation  is  used  for  P^  fitted  separately  to  the 

saturated  vapor  pressure  data  of  Cragoe^  for  the  range  from  the  triple  point 

7  8 

to  373K  and  of  the  mean  of  the  data  of  Beattie  and  Lawrence     and  Keyes 
above  373K.     The  equation  so  obtained  is 

loe    P  =    -  +  B  +  CT  +  DT2  +  ET 
•LOge    s  T 

T  in  kelvin  (IPTS68),  P  in  atm. 

-  A  -  -3684.7798 

-  B  =  20.428787  (10b) 

-  C  -  -.02893289 

-  D  -  3.4798128xl0~5 

-  E  -  -9.2219845x10' 9  . 

Equation  (10b)  and  eq     (3)  define  the  coexisting  phases  for  this  report  and 

all  properties  over  the  temperature-pressure  range  of  the  thermodynamic 

surface,  subject  to  the  low  temperature  boundary  of  the  melting  solid. 

The  thermodynamic  surface  is  consistant  with  the  following  values  for  the 

parameters  at  the  triple  point, 

T    =  195.48  K 
t 

P    =  .06063  bar 
t 

gt  =  .00006382  g/c3,        =   .73374  g/c3, 


P 

and  at  the  critical  point, 


T    =  405.4  K 

c 

P    =    113.04  bar 

c 

Pc  =  .2350  g/c3 


7a 

The  relationship  between  the  pressure  and  temperature  of  the  melting 
solid  was  calculated  by  means  of  the  Clapyron  equation,  using  the  latent  heat 

of  fusion  and  the  specific  volumes  for  the  saturated  liquid  and  solid. 

9 

For  the  latent  heat,  the  value  reported  by  Overstreet  and  Giauque    was  used;  for 
the  specific  volume  for  the  solid  at  the  normal  melting  point,  the 

value  reported  by  McKelvey  and  Taylor"^  was  used,  and  for  the  corresponding  specific 
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11 

volume  of  the  liquid     the  value  reported  by  Cragoe  and  Harper      was  used. 
The  Clapeyron  equation  is  the  relation 


dp,  an 


where  the  quantities  u'   and  u  refer  to  the  specific  volumes  of  the  liquid 
and  solid,   respectively,  and  L  refers  to  the  latent  heat  of  fusion.  From 
the  above  data,  the  quantity  — r- — —    ■-     4x10  ^  atm  ^ ;  and  Eq  (11)  can  be 
integrated  to  yield 


T  =  Tg  exp[4xl0  5 (P  -  Pg)],  (12) 


where  T    and  P     are  the  triple  point  values.     (The  differences  between  the 
s  s 

triple  point  values  and  those  of  the  normal  melting  point  are  negligible.) 
Also,  since  P^  at  the  triple  point  is  very  small,   the  relationship  can  be 
simplified  to 


T  -  195.48  exp  {AxlO-5  P  (atm)}.  (13) 

J.       Computer  programs. 

The  properties  calculated  from  the  equations  presented  in  this  report 
have  been  compared  with  the  various  thermodynamic  measurements  reported  for 
ammonia  in  the  comprehensive  review  in  reference  1.     In  almost  all  cases  the 
agreement  is  within  the  experimental  accuracy  of  the  data.  Comparisons 
of  results  with  P,p,T  data  are  given  here  in  figures  2  and  3.     It  was 
established  in  reference  1  that  for  most  of  the  vapor  phase  and  for  the 
coexisting  phases,   the  calculated  enthalpies  are  accurate  to  within  0.1%. 


To  facilitate  the  application  of  these  results  we  present  in  Appendices 
I  and  II  computer  programs  with  which  the  various  properties  of  ammonia  can 
be  calculated.     Two  such  programs  are  included: 

Computer  program  I  in  Appendix  I  refers  only  to  certain 
properties  for  the  coexisting  phases,   including  the  latent 
heat  of  vaporization,   the  vapor  pressure  of  the  saturated 
liquid,   the  densities  of  the  saturated  vapor  and  the 
saturated  liquid.     The  dependent  variable  for  these  is 
either  the  temperature  or  the  corresponding  pressure  of 
the  saturated  liquid. 

Computer  program  II  in  Appendix  II  is  a  general  program  for 
all  the  thermodynamic  properties  discussed  in  section  2  of 
this  report.     The  dependent  variables  are  either  pressure 
and  temperature  or  density  and  temperature. 
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Figure  Caption 

Fig.   1.       P-T  Schematic  of  the  data  included  in  this  correlation.  The 

various  polygons  represent  the  range  of  the  individual  sets.  The 
numbers  refer  to  references  in  this  report  for  the  sources  of  the 
various  data  sets. 

Figs.   2,3.   Comparisons  of  calculated  values  with  P,o,T  data.     The  results  are 
presented  as  fractional  differences  in  percent  vs  pressure  in 
atmospheres.     The  numbers  to  the  left  of  the  various  plots  refer  to 
data  references  in  the  text.     Figure  2  refers  to  pressure  deviations 
for  densities  less  than  1.5  times  the  critical  density.     Figure  3 
refers  to  density  deviations  for  values  of  density  greater  than  1.5 
times  the  critical  density. 
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Appendix  I 

I.J 

Listing  oft  Fortran  program  I: 

Temperatures,  Pressures  and  Densities  along  Saturation  Curve. 

A.  The  saturation  pressure  as  a  function  of  temperature  and  the 
saturation  temperature  as  a  function  of  pressure;  where  T  is 
in  degrees  K,  and  p  in  atmospheres. 

Function  PS(T) 
Function  TS(P) 
Function  DPSDT(T) 

B.  The  saturation  liquid  and  vapor  densities  as  a  function  of 

3 

Temperature-:    (densities  in  gram/cm  )  . 

1.  SUBROUTINE  DS(T,DL,DG)   to  calculate  approximate  densities 
from  a  simple  polynomial  to  be  used  as  initial  guesses  for 
the  subroutine  DFIND  which  will  improve  upon  these  by  an 
iterative  process. 

2.  SUBROUTINE  DFIND (DOUT ,P ,DGUESS ,T ,DPDD)  where  P,T,DGUESS 
are  inputs,  DOUT  and  DPDD  (8P/8p)  are  outputs  and  which 
also  uses 

3.  Function  PRES(D,T)  which  calculates  the  pressure  in  atmospheres 
at  any  point  D  and  T 

SUBROUTINE  QQ(Q,L,K,T ,D)  which  calculates  "Q"  or  the 
8L+KQ 

derivatives  — at  an^  P°int  T  and  D. 

3p  9t 

C.  SUBROUTINE  HV(HH,T,DL,DG)  which  calculates  the  heat  of  vaporization 
as  a  function  of  T  and  the  previously  calculated  liquid  and  vapor 
densities.     The  HV  calculated  here  will  be  in  dimensionless  units, 
and  to  obtain  results  in  joules  per  gram  multiply  by  .4882  T  or  in 
BTU/lb.  multiply  by  .210027  T  (T  in  both  cases  in  kelvins) . 


D.       Units  and  conversion  factors:  a  package  of  subroutines  is  included 

which  will  convert  Temp,   in  °F  to    K  or    K  to  °F,  and  P  in  PSIA  to 

3  3  3 

atm.   or  atm.   to  PSIA,  and  densities  in  lb/ft     to  gm/c     or  gm/c 

3  1  '  ' 

to  lb/ft  . 


1.3 

FUNCTION  PC(T) 

PL--3684i  7798/T  +  20,+2  8787  +  T*(  -  .n2  893  289+T*'(  3,4798l28r-5-T*9.2  2l984 

1  ^cr-O  )  ) 

PS=FXP(PL) 

RETURN  v 
END 


FUNCTION  TS(P) 
0=1  .  VP 

T  =  252. 10703   +.  P* ( 5 ,66586  36+P* ( -. 12457077+P*( .OOl3  5987-P*,OOnoo5319 

142)))    +  Q*(-i  ^.37P67R+Q*(  1  .^2^1  5-Q*. 04906)  ) 
nTr>p  =  i  .  /nPSPT  (  T  ) 

5   PX  =  P  ^ ( T ) 

TF(  ARS(  l.-P/PX)  .LT.l  m*-5)  TO   ].  0 

PT=  (  P-Dy  )  tfnTDP 
T=T+PT 
00   TO  * 
to  TS=T 
Q  c  T 11 R  N 
F  M  0 


FUNCTION  OPSDT(T) 

DF=-.0  289  3  28  9+3684.779  8/T/T+6.9  596256r-5*T-2.766  5  9535F.-8*T*T 

r>PSOT  =  nF'*P  s  (,T) 

RETURN 

FNO 


SURROUTTNF  0S(T>PL,D0) 
OOUPLF   PRECISION  T,OL»nr-, 
n  c  =  .  ?  0 : 5 

TF( T-4n^.4)  in,fln,on 
10   nTri  ,-T/4nt;.4 

T  c ( T , L  T . 1 n  ^ .  )    GO    TO  iOO 
T  V=DT**.^5 

T2=DT**.541 
T-5  =  T2*T2 
T4-nT^2.o^c,^ 

nL  =  T4  +  Tl^(',.i  i  ^  —  1  ,Ano7«T?-,Roflo?*T^) 
nr.  =  T4-T  V*  (  2  .  11  7  +  1  .  ]  ?«>0*T?  +  .  ^7?5**T**  ) 
DI_  =  0C* ( DL+1 . ) 
0  g  =  n  r  *  (  P  r-.  +  i  .  ) 
RFTURN 

po  0L-nr 

D^=.DC 
P  r  T 1 J  P  M 
no  DL=0. 
on  =  0 . 

R  F T'  IRN 
1  On      n  =  OT 

nL=. 387 13 1 -.00096947 /Q+Q* ( 1 .15 138  7 5+Q*< -1.494  3 106+U*1.1 18  3325) ) 
DG=. 08 28 86  7+. 0009 5867/Q+Q* ( -.6  0039  5  34+0* ( 1 .443  459  4 
1      -0*1  .  167*60R  )  ) 
RETURN 
END 


1.4 

SURROUT iNr  DFT  Nn ( ^OUT  »P  »n , T , nPn ) 
COMMON   / UQQQ /  00 »Q 1 »Q2 » 01 0 »Q20 »Q1 ] 

D  0  =  D 

I F ( T .LT.223 .  1    •  A  N  H  .    n.r-T.,01    .  ANO.    D. LT..7)    ^0=.7  3 
1  =n 
o  L=L+1 

I  F  (  no.LF.O.  )    nO=l  .  F-R 
I  F  (  f)  D  .  0  T  .  .  o  )    p  n  =  .  o 

TALL  QQ(QO»0»OtT»r>P) 
CALL  QO  (  Ql.  ,  0»  1  9  T  »0D  ) 
CALL  00 (0?  »n  ,? ,  T ,DD) 
npn=np^p (nn,Ti 

T  F  (  npnx  .  LT  . «;  0.  )    PP^X  =  r^. 
PP=PPFS ( OD , T ) 

TF(  ARS(  ]  ,-PP/P)  .LT.l.F-^)    r.O    TO  ?Q 

X= ( P-PP ) /PPPX 

OD=OD+X 

I F ( L , LF«  ?0 )    GO   TO  o 

°0  CONTINUE' 

R  F  T  U  f?  N 


FUNOT I  ON   T K  ( TF ) 

R  ^  T  I  !  P  N 
FMn 

F  U  N  C  T  T  0  M  P  A  T  M  (  P  S  I  A  ) 
P  A  T  m  =  p  $  I  A  /  1  4  .  *  ^  * 
P  F  T 1 1  R  N 
FMD 

FUNCTION   P    r  r  (  n  L  P  r  r  , 

nrn r  r  =  n  l d  r  r  /  a  ?  .  4  ?  R 

RETURN 
FNH 

FUNCTION   DLBCF  (  DCCr  ) 
n  |_  p,  r  f  -  p  r,  r  r  *  a  ?  .  4  ?  R 
PFTURN 
FMD 

FUNCTION  PSTA(PATM) 
P  $  T  A  -  P  A  T  M  *  1  4  .  A  o  6 

PFTI  IRN 
rMD 

FUNCTION    T  r ( T  K  ) 
TF=l  .  P  *  T  K  -  4  r^  .  ^>  7 

PFTURN 
F  ,M  P 


SURROUT!Nc   00 { 0 , L , K , T , l  ) 

n  T  M^NS  t  ON   P  (  5  »  5  )  »        1  0  )  >  TT  (  7  )  ♦  A  (  44  )  »  T  T  (  44  )  ♦  JJ  (  44  )  ?  Al  f  22  )  »  A2  (  22  ) 

foutVal-N^1"    (  &1  (  1  )  *  a  (  1  )  )  ,  (  a?  (  i  )  ,  a  (  ?i  )  )  1 

nATfl    f  t/6*1 ,6*2  » 6*3 » 6*4  *4*5 ,4*6 ,4*7 ,4*8 ,4*9/ »N/44/ 

^/\Ta    JJ/i»2,  3  ,4,5, 6»1»2»3»4,  5  ,6,1, 2,3,4, 3,6»l»^>2»,:+,i,6>l»2»3»/-+»l» 
1   2  ,  3  ,  4  »  1  »  2  »  3  *  4  »  1  »  2  ,  3  *  4 , 1  ,  2  ,  3  »  4  / 

DATA  P/2.>2.»3.»4.,5.»2.,4.,6.»12.,20.»3.»6.,12.»24.,6o..,4,,l2., 
1    24.  ,48.  » 120.  ,5.  ,20.*, 60.  ,  l2o.  ,240./ 

DATA   Al /-6. 469043-56    ,-13.29562588    ,-8.121177092    ,-6.969004355  , 
]    -9.736580235  .  ,3.481664262    ,8.81  0044576    , -5  .  0789  56.87  1  ,-68.2615834 
22    ,-74.  /271569b    ,49.  75165418    , -14.^6  71  5637    ,-  1  o  .  u  6  7  ^ri  2  86  ,361,f>19^ 
3765    ,1327.827022    ,  1484.264330    , -d 2. 22 9 12 2 94    ,20.17085672  , 
4   75.04957400    ,-2103*945194    , -7b 76 . 1 00794    ,-8334.874642  / 

DATA   A.2/43. 99847596    ,-9.2  7  733  7672    ,-409.0296415    ,6212.282252  , 
122341.80033    ,23618.79174    ,1072.479955    ,-10816.10642    ,-38259. 3441 1 
2,-38233.53  40  0    ,-1471.40.1315    ,11195.13872    ,38544.62819  ,35887.29465 
3    ,1046.234130   ,-6365.756670    ,-2131^.81^31    ,  - 1 8  1  62  .  o  949  7    ,~30 5.3008 
41  17    ,1  532.061605    ,502 T.6'96209    ,^812.^6^15?  / 

0^0. 

T  P ( L  +  K )    1  ?  ,  1 4  ,  1 p 
1 ?  SATURN 
1A  U=SO0./T 

r-u— 1  .?^^?-4o77Pr;0 

fF(DARS(r).LT.l.n-q.)    r  =  l  ,  n  -  8 

TT ( 1  )  =  1  . 

00    IS  1=7,7 

lis   TT ( I ) =TT ( T-i )*C 

T  r  (  n  .  |_T  .  1  .  F-8  )  0=1.^-3 
DD( 1 ) =1 . 

.no  16  T  =2 » i 0 
i 6  nn ( t ) =dd ( t -i ) *r> 
i fl   DO   200  M=i ,N 

t  =  T  I  (  M  ) 

J=JJ (M ) 

T  F ( J-L-1  )    ?oo, ?0 
?  °    I  F  (  L  )    ?00 , ?4  ,  ?8 
?4   QT  = l  . 

oo  to  so 
?  p   0  T  =  P  (  L  ,  L  )  /  ?  . 

r-n   TO  50 
on    jr(D    -7  00,"*  ?,40 
^?   OT=TT ( J ) 

<^0    TO    ^  O 

6.0   0T  =  P  (  L  ,  J-i  }  *TT  (  J-L) 
so    IF(K)    2  0 "~  ,  A  o  ,  7  p 
6^   QT=QT*on ( ] ) 

on   TO   1 0  0 
70   P=i . 

no   7?   vsm  =  1  ,  k' 

-79    D  _!}-«-  (  J  -MM  ) 

O  T  =  0  T  *  P  *  0  D  (  t  -  y  ) 

1nO    Q-O+ft  (  V.  )*QJ 

?nnroNTTMI,F 
RrTUr?N 
FNO 


r  i  0 r-J     '-'V  (  T  ,  r,|  ,  nr.  ) 

C  a  LL  00 ir«':  ,  Q,Q,  T  ,nL  ) 

r  ALL  00  (  01  »  ^  »  l  ♦  T  *>!..  l 

TALL  00  (  0?  » i  ?.n  ,  T  »  0L  ) 

r  A  LL  0 0  (  S r  »  ^  » n  .»  T  »  n r-  ) 

r  ALL  00  (  Si  ,  r, ,  i  ,  T  »  ^ ) 

CALL  OQ ( S  ?  » 1  »0  »  T  »  DO ) 

HV  =  U*  (  ^r,*5  ?-nLi;C?  )    +  (  SO  +  nr,'Sl  )    ~   nL*  (  UO  +  ^L*vj  ]  ) 


A   SAtfPL^   PROGRAM   To  r  A  l  r  U  L  A  T    »    coR    r-NSTAtVrc,    Tnr  b  ATURa  T  I  OiNi   TF/.P^R  a  T.UR 
IN  Dc0  p»    DENSITIES    iN  L  5  /  CU   F  T    ANp  HP  AT   OF   VAPOR  t  Z  aT  I  ON    rN   pTU/L?   aT  a 
r,yV="N   PRrSSURr    T  N   PSta   rO'jLn  AS   FOLLOWS  - 

P'FA?    i  ,  £  S  ! 
i    f  0  p  m  a  T  (  *  P  t  C  .  ^  ) 

D_D  £  TU  {  P5  I  } 

T  =  T  S  (  P  ) 

r  A  l  L   D  S  (  T  »  X  0  L  »  X  ^  ^  ) 

CALL  nPTMO(nL,P,XPL»T.nX) 

'"ALL        I  Nn  (  nr, ,  p  ,  xnr. ,  T  ,  nY  ) 

uu  =  wv_(  T»      » fS'G  )*.•'*>  i  30?  7*  T 
nL  =  OLa^F(n[_) 

TrTF(T) 

PRINT    T»   P5I »TT ,DL ,D^,HH 
c  tod 

F  N  ^ 


Appendix  II 

Listing  of  Fortran  Program  II: 

Routines  for  calculating  at  any  point  on  the  thermodynamic  surface. 

Subroutine  FZ (T , F,DF , DDF , DTF)  which  calculates  for  a  given  temperature 

T  in  deg  K,  the  ideal  gas  functions  for  the  free  energy  F,  entropy  DF, 

C    DDF,     and  internal  energy  DTF,  all  in  dimensionless  units. 
P 

Functions  for  calculating  internal  energy,  enthalpy,  entropy, 

CL  and  C    as  functions  of  D,T  and  ideal  gas  functions,  all  in 

P  v  ' 

dimensionless  units. 

Phase  finder:  SUBROUTINE  PHASE  (P , T , DL  ,DG ,K)  which  will,   for  any 

point  P  and  T  supplied,  return  a  value  of  K  as  follows: 

K=0  point  is  on  sat.  curve,   guesses  for  DL  and  DG  returned 
K=l  T<T  ,  vapor  phase,   guess  for  DG  returned,  DL  set  =  DG 
K=2  T<Tc,  liquid  phase,   guess  for  DL  returned,  DG  set  =  DL 

K=3  T>T    DG=DL ,  guess  returned. 

c 

K=-l  error  return  for  T  <  or  P  <  0.   or  for  either  P  missing 

and  T  >  T  . 

c 

for  T<T  ,  and  P  missing  or  =  0. ,  P        will  be  calculated,  and  the 

C  S  3 1 

appropriate  densities  found.     Similarly,  for  P<Pc  and  T  missing  or  =  0., 

T        will  be  calculated, 
sat 


II. 2 

SUPROUTl NF   FZ( T ♦ F » OF , DHF > DT F ) 
OOUPLF   PRECISION  T 
D I MFNS I ON  A( 1 ?) 

DATA    A/-3. 872727, .64463/ 24, J. ^3871;  9, -,nn21  3  76925, .ob89nbJ3F-I>. 

1  -  .  2  4  0  8  f>  I 4  9  F  -  7  ,  .  3  6  8  9  3  1  7  5  F  - 1 0  »  -  •  3  !>0  3  4  6  6  4  F  - 1  3" »  .  2  0  5  6  3  0  2  7  F  - 16, 

2  -.68  5  342E-2  0,.?9  3  9242  7E-24,0./ 
U  =  T 

I F ( T.LF.l 0. )  U=10. 

F    =  a  (  1  )  *  A  LOr,  (  U  )  +  A  (  ?  )  /U+A  (  ~  )  +  A  (  a  )  *U:+  A  (  <=>  )  *U*U4-A  f  6  )  *U*  *^  +  A  (  7  )  *U*#4 
1    +A  (  8  )  *U,#*'5  +  A  <  0  )  *U  v  *6+A  MO)  *U**7.+  A  (11)  *U**P 

HF  =  A  (  1  )  *  (  ALOG  (  U  )  +1  .  )    +   AO)    +   2.*A(4)*U  +    3  .  *  A  (  5  )  *'U  *'U   +  4.*A(6)*U 
+    5*^A,(  7  )*il**4    +    6  •  *A  (  °  )  *U*  *r-    +   7  .  *  a  (  o  )  +    p  . -""A  (  1  ^  )  *M;; --'7 

?    +    Q  ,*A (  1  1  ) *U**8 

DnF  =  A  (  1  )  /U  +  2  .*  A  (  A  )  +6  .  *  A  (  5  )  *  U+  1  2  .  *  A  (  6  )  *U*U4  2  0.  *a  f  7  j  lfU**?  +  3  0  ,*A  (  fl  )  * 
1    U**4   +   42  .  *.A  C  9  )  *U**5   +   56.*A (  1  0  )  *U**6   4   7  ?  .  *  A  (  1  1  )  *U**7 

DDf=-ODF^U 

D  F  =  -  D  F 

DTF  =  nF  +  F-l  . 

RETURN 

END 

subROuti nf  Phase ( p * t ».d.l , dg » ^ ) 
i  F ( P  )  10,10,20 

iO    if(T-i°-.  )  on,io,T9 
17    T  F  (  T  -  A  0  r  .  5  )    1  A  ,  o  0  ,  o  0 
1  4   P  =  P  5  (  T  ) 
16   CALL   P  5  (  T  »  n  L  »  n  O  ) 

K  =  0 

RETURN 
20    I  F  (  T  )    7.2*  22  »  3,0 
22    I  F  (  P-l  1  1  .  5^  )    24*  ?4  »"o0 
?  4    I  F  (  P  -  •  0  A  }    9  0  >  ?  6  »  ?  6 
?  6   T  =  T  S  (  P  ) 

GO   TO  16 
30    TF(T-4  0r>  ^?»<0>^6 
3?   PX-PS (  T ) 
34  CALL   D S ( T ♦ D L  » D G ) 

IF(P-PX)    36  »•  14»40 
?6  DG=DG*P/PX 

DL  =  DG 

K  =  l 

RETURN 
4  0  DG=DL 
K  =  2 

RETURN 

50   D  =  P*(..0002107-5.19F-8*P-2.289F-7^T  ) -1  .2823  F-6*T 
D  =  D  +  P / T/4 • 8 1 8 
I E ( D • LT •  1 »F-5 )    D= .00001 

PL  =  D 
DG  =  DL 
K  =  3 

RETURN 
QO  DL=P. 
DG  =  0. 
K  =  -l 

R  FTURN 
END 


II.  3 


SUPROUT  T N^   THERMF ( D » T » rP , C V  « F , H » S ) 
COMMON   //'JQQQ/  Q 0 0  » 0 Q  1  » Q  0  2  » Q 1 0  » 0  2  0  » Q 1 1 
U  =  «500  ./T 

CALL  QQO(QOQ,n»0»T»D) 
CALL  QQQ(QO] >0»1 »T,D) 
CALL  QQQ(QO?»0»2»T,D) 
CALL  QOQ ( Q 1 0  » 1 »  0  »  T  »  D ) 
CALL  QOU (Q?0»2*0,T,D) 
CALL  QQQ (Ql 1 » 1 ♦ 1 9 T ,D ) 
CALL  FZ(T,FT»SIG»CPIG»FIG) 
ENFR  =  D*U  »-  Q 1  0  +  FIG 

FNITH='FNER   +    1.    +   D#OC)r>   +  D*D*QOl 

FNTR  =  S  T  G  +  D*U*Ql  0-0*0  0  0- A  LOG  (4,81  R*D*T) 

CV=CPIG-1 .-D*U*U*Q?0 

rP=  ( i  ,  +  n*  ( ooo-u*ui  o+n*  (  ^o  i-'J*oi  i )  )  j  **2/  ( i  ,  +  n*  (  2 .  *ubo+r>*  (  4  ,*uo  i  +  d* 

1   Q02 )  )  )  +CV 
RETURN 
END 


thp  following  samplf  program  will  calculate  the  density  and  fivf  thermo- 
dynamic FUNCTIONS  AT  A  POINT  (P,T)  5UPPLIFD.  P  SUPP^-IFD  I N  PAR,  T  I N  DFG  C 
h  F  N  S  I  T  Y  IN  GM/CC*  rP»CV»S  IN  JUIJL  FS  /  Givi  DFG  r*  AND  iNIfRNaL  ENERGY*  ENTHALPY 
IN  JOULE S/GM. 

1  READ   2  »    PB AR  »  T C 

2  FORMAT  (  2F10.3  ,F10..6  »2F10.3  ,3F10.5  ) 
P'=P.RAR/]  .0]  ^2^ 

IFtP.Lr.C.)  STOP 
T  =  Tr  +  ?7%  1  S 

CALL  PHASF(P»T»DL»DG»K) 
I F ( K.FQ.-l )    GO   TO  1 
CALL   DF TND ( D »P » DL , T » DX ) 
GALL   THFRMF( D»T  »CP  »CV,F»H»S ) 

tP=CP*«48R? 
CV=CV*. 4RR2 
S  =  S  *  •  4  R  R  ? 
H=H*.4RR?*T 
E=F*.4RR?*T 

PRINT  ?>PBAR>TC»DvH>F»S»cP»CV 

GO   TO  ] 

END 
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